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Algorithm of ensemble pulsar time 
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Abstract An algorithm of the ensemble pulsar time based on the Wiener 
filtration method has been constructed. This algorithm has allowed the sep- 
aration of the contributions of an atomic clock and a pulsar itself to the 
post-fit pulsar timing residuals. The method has been applied to the tim- 
ing data of the millisecond pulsars PSR B1855+09 and PSR B1937+21 
and allowed to filter out the atomic scale component from the pulsar phase 
variations. Direct comparison of the terrestrial time TT(BIPM96) and the 
ensemble pulsar time PT cns has displayed that the difference TT(BIPM96) 

- PT ens is within ±0.4 /is range. A new limit of gravitational wave back- 
ground based on the difference TT(BIPM96) - PT ens was established to be 

n g h 2 ~ io- 10 . 
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1 INTRODUCTION 

The discovery of pulsars in 1967 ( Hewis h et al. 1968j) and millisecond pulsars in 1982 
IjBacker et al. 198 21 and consequent observations had shown clearly that their rotational 
stability allowed them to be used as astronomical clocks. 

In this paper the author presents a method of forming the ensemble pulsar time scale 
(PT). The method is based on the optimal Wiener filter. In Sect. principles of pulsar 
timing are described with regard to time scales. Sect. contains a theoretical algorithm 
of the Wiener filter and construction of the ensemble pulsar time scale. Sect. 0] presents 
results, Sect. El discusses an application of the algorithm to timing data of pulsars PSR 
B1855+09 and PSR B1937+21 ( |Kaspi, Taylor fc Ryba 19941 ). 
* E-mail: rodinOprao . psn . ru 
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2 PULSAR TIMING 



An observer which is situated on the Earth rotating around its axis and moving around 
the Sun receives with a radio telescope a pulsar signal during an integration time to 
obtain sufficient signal-to-noise ratio. Time of arrival (TOA) of the integrated pulses 
are measured with the observatory freguency standard (e.g. H-maser) by the maximum 
of the cross-correlation between the integrated pulse and the pulse template. The ob- 
tained topocentric TOAs tn are in the scale of the local frequency standard and there- 
fore required to be converted to the barycentric time scale via the following expressions 
dManchester fc Taylor 1977| |Doroshenko fc Kopeikin 1990| ): 



where Ati is the correction of the local scale to the universal coordinated time UTC. 
The international atomic time (TAI) differs from UTC by k integer number of seconds 
introduced to reconcile the lengths of day measured by an atomic clock and the Earth 
rotation. TAI is related with the terrestrial time TT by the eq. (JIJ, where the constant 
shift 32.184 sec prevents a jump between ephemeris and atomic time. Since a second 
of TT has various lengths depending on the position and velocity of the Earth in its 
orbit then a transformation from TT to TB scale is required which is performed on 
the basis of the paper by Fairhead & Bretagnon (1990). Once converted from TT to 
TB the topocentric TOAs need to be reduced to the barycentre of the Solar system 
(SSB) according to the following transformation formula ( [Manchester &: Taylor 19771 
|Doroshenko fc Kopeikin T 990 ) : 

T = t-t + A R (a, S, [i a , [is, tt) - A_r(x, e, P b , T , oj) - D/f 2 + A rol + At c i ock , (2) 

where to is the reference epoch, t is the pulsar topocentric TOA in TB scale, T is TOA 
at the Solar system barycentre in TB scale, A#(a, 5, \x a , /is, tt) is Roemer delay along the 
Earth orbit, Ar(x, e, Pb, To, lo) is Roemer delay along the pulsar orbit, D/f 2 is dispersive 
delay for propagation at the frequency / (corrected for the Doppler shift) through the 
interstellar medium, A re i is time corrections due to relativistic effects in the Solar system 
and the pulsar orbit, Ai c i oc k is the offset between the observatory frequency standard 
and the terrestrial time. 

Time of arrivals at the SSB are then used for calculation of the pulsar rotational 
phase (in cycles) 



where Nq is the initial phase at epoch T = 0, v, i> are the pulsar spin frequency and 
its derivative respectively at epoch T — 0, s(T) is phase variations (timing noise). The 
fitting procedure includes adjustment of the parameters No, v, i>, a, S and so on, to 
minimise the weighted sum of squared differences between N(T) and the nearest integer. 



UTC = t n + An, TAI = UTC + k, TT = TAI + 32.184 s, 



(1) 



N(T) = No + vT+ -vT 2 + e{T), 



(3) 
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Usually the pulsar rotational phase residuals are expressed in units of time St = 5N /v. 
In this paper we deal with the only part of the residuals that includes the variations of 
the clock offset At c i oc k(T) . 

When comparing different realisations of atomic time scales between each other one 
can see that they are dominated by flicker frequency noise on intervals of a few months 
and random walk in frequency on intervals of years IjOuinot 1988JI . i.e. in the frequency 
domain the clock variations have power spectrum of form l/u> n , in the time domain clock 
variations can be expressed in the polynomial form 

A* c i ock (T) = c + cT + hr 2 + \cT z + . . . . (4) 
2 o 

One can see that the appearance Ai c i oc k in eq.JSJ results in a coupling between pulsar 
and clock parameters: 

N(T) = ^ + (1 + c)/T + ~(/c + (1 + cff)T\ (5) 

where T is the ideal time scale, /, / are the pulsar frequency and its derivative not 
subjected to the influence of the clock parameters. For this reason one should use TO As 
expressed in the best available time scale TT (jGuinot fc Petit 199 1|> . 

3 FILTERING TECHNIQUE 

Let us consider n measurements of a random value r = (n, r%, ■ ■ ■ , r n ) are given, r is a 
sum of two uncorrelated values r = s + e, where s is a random signal to be estimated 
and associated with the clock contribution, e is random errors associated with the fluc- 
tuations of pulsar rotation. Both values s and e should be related to the ideal time scale 
since a pulsar on the sky "does not know" about the time scales used for their timing. 
The problem of the Wiener filtration is concluded in estimation of the signal s if the 
measurements r and the covariances (si,Sj) and (ei,Sj), — 1,2, ... ,n) are given 
IjOubanov 19971 [Vaseghi 2000| ). The optimal estimation of the signal § is expressed by 
the formula l|Gubanov 19971 |Vaseghi 2000| ) 

§ = Q sr Q^r = Q^Q^r = Q SS {Q SS + Q ee ) _1 r, (6) 

where Q rr , Q ss are the covariance matrices of the noise data r and signal s respectively, 
Qsr, Qrs are the cross-covariance matrices between r and s. The covariance matrix Q ss is 
calculated as cross-covariances ( h ri, l Tj) — (si, Sj) , (fc, Z = 1,2,..., M; i,j = 1,2, ... , n), 
M is a total number of pulsars. In the formula © the matrix Q^r serves as the whitening 
filter. The matrix Q ss forms the signal from the whitened data. 
The ensemble signal (time scale) is expressed as following 

M(M-l) 

2 2 M 

§cns = M(M 1) £ fc Q--£^Q--*r, (7) 

v ; k=i i=i 

where Wi is the relative weight of the ith pulsar, Wi oc 1/ci, is the root-mean-square 
of the expression 'Q" 1 • l r. 
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Fig. 1 The barycentric timing residuals of pulsars PSR B1855+09 (dots) and 
PSR B1937+21 (continuous line). 



4 RESULTS 

The method described in the previous section has been applied to the pulsar timing data 
of PSR B1855+09 and B1937+21 flKaspi, Taylor fc Ryba 19941 ). Though these data are 
regular they are unevenly spaced, therefore a cubic spline approximation has been applied 
to make them uniform with sampling interval 10 days. Such a procedure perturbs a high- 
frequency component of the data and leaves unchanged a low-frequency component which 
is of interest. 

The common part of the residuals for both pulsars (251 TOAs) has been taken within 
the interval MJD = 46450 -j- 48950. Since the residuals after the procedure of dropping 
their ends have the different mean and the slope they have been quadratically refitted 
for consistency with the classical timing fit. The residuals after all treatments described 
above are shown in Fig. 1. 

According to ( |Kaspi, Taylor k Ryba 1994| ) the timing data of PSRs B1855+09 and 
B1937+21 are in UTC time scale, hence the signal to be estimated is the difference 
UTC - PT. Fig. 2 shows the signal estimates based on the residual TOAs of pulsars 
PSR B1855+09 and PSR B1937+21. The ensemble signal and the difference UTC - TT 
display similar behaviour. 

5 DISCUSSION 

For calculation of the fractional instability of a pulsar as a clock a statistic a z has been 
proposed (Tay lor 1991| ). A detailed numerical algorithm for calculation of a z has been 
described in the paper ( Matsakis, Taylor, Eubanks 1997| ). 

Fig. 3 presents the fractional instability of PSR B1855+09, PSR B1937+21 and TT 
- PT ens . The theoretical lines of a z | |Kaspi, Taylor fc Ryba 1994| ) in the case when the 
timing residuals are dominated by the gravitational wave background with flgb 2 = 10~ 9 
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Fig. 2 Combined clock variations of UTC - PT in interval M JD = 46450 
48950 estimated using the optimal filtering method based on the timing residuals 
of pulsars PSR B1855+09 (thin line), PSR B1937+21 (dashed line), ensemble 
UTC - PT cns (dot-dashed line) and UTC - TT (solid line). 



and 10 10 are plotted in the lower right hand corner. One can see that a z of TT - PT ens 
crosses line fl g h 2 = 10 -9 and approaches flgh 2 = 10 -10 . 

The fractional instability of the TT relative to PT ons is at level of 10 -15 at 7 years 
interval and almost one order better than the fractional instability of the pulsars PSR 
B1855+09 and PSR B1937+21. It is expected that reliability of TT - PT ens estimation 
will grow up by increasing the number of pulsars participating in PT ons as M(M — 
l)/2 (the number of cross-correlations). Currently the accuracy of the filtering method 
described above without contribution of the uncertainty of TT algorithm is estimated at 
level 170 ns. This uncertainty is obtained as root mean square value of the data points 
taken within the smoothing interval of span m. The span m was calculated from the 
equivalent bandwidth of the low-pass filter applied to the ensemble data for more easy 
comparison with UTC - TT line. The uncertainty of PT cns may, in principle, reach the 
nanosecond level if to use all the observed highly-stable millisecond pulsars. 

The method proposed can not distinguish the 2nd order polynomial trends in the 
reference clock and the pulsar phase due to pulsar period slowing-down. However, this 
is not a problem if to consider the timing data at more long intervals and process them 
off-line. Under such processing the long-term details appears as the data span is increased. 

The low fractional accuracy of P mentioned in the paper IjGuinot fc Petit 199 1|) pro- 
duces no disadvantages when processing off-line since no prediction of the pulsar rota- 
tional phase is performed. However if one does need to predict a behaviour of the concrete 
atomic scale variations, e.g. UTC, then this can be done on the basis of the UTC - PT cns 
data by using standard forecasting methods for the time series, e.g. the auto-regression 
method with reservation that only relatively short-term variations without quadratic 
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Fig. 3 The fractional instability a z for pulsars PSR B1855+09 (solid line), PSR 
B1937+21 (dashed line) and TT - PT cns (dot-dashed line). 

trend are forecasted. Under such an approach the unsatisfactory fractional accuracy of 
the spin period derivative does not play significant role since the phase variations are 
predicted rather than an absolute value. 

The proposed filtering method can be applied in "inverse" form: one pulsar and a few 
reference clocks. In such case it is also possible to separate the pulsar timing noise and 
the clock variations relative to the ideal time scale rather than to obtain a simple clock 
difference. 
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